A Stochastic Liouvillian Algorithm to Simulate Dissipative Quantum Dynamics 

With Arbitrary Precision 
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An exact and efficient new method to simulate dynamics in dissipative quantum systems is presented. A 
stochastic Liouville equation, deduced from Feynman and Vernon's path-integral expression of the reduced 
density matrix, is used to describe the exact dynamics at any dissipative strength and for arbitrarily low tem- 
peratures. The utility of the method is demonstrated by applications to a damped harmonic oscillator and a 
double-well system immersed in an Ohmic bath at low temperatures. 



Relaxation dynamics and fluctuations of a quantum system 
in contact with a thermal bath are improtant for many different 
areas of chemistry and physics, such as reaction-rate theory 

ultrafast phenomena in chemistry and biochemistry 
tunneling at defects in solids [[|] and quantum optics. While 
some limiting cases permit Markovian or weak-coupling ap- 
proximations to the dissipation effects, the majority of others 
do not. Problems involving coherence [|[] are prime examples 
in which the Markovian approximation fails. 

The exact treatment of dissipative dynamics presented here 
is based on the path integral influence functional formalism 
developed by Feynman and Vernon for interacting sys- 
tems and extensively used in dissipative quantum mechanics 
The main difficulty in evaluating the complex-valued 
path integral involved is that the integrand is not a local func- 
tional of the paths. Except for the case where the memory 
effects decay on a short time scale deterministic integra- 
tion methods fail, making it prohibitively expensive in both 
computation time and memory requirements. To date, the 
only exact numerical approach to dissipative quantum dynam- 
ics with arbitrary memory time has been the dynamical quan- 
tum Monte Carlo (QMC) method. This method, however, is 
severely limited by the sign problem, i.e. the signal-to-noise 
ratio of the numerical result tends to zero exponentially with 
the time scale of the problem. 

In this Letter, we deal with the problem with non-local 
influence functionals arising from arbitrary long bath mem- 
ory. We recover fully time-local dynamics for the ubiqui- 
tous case of Ohmic dissipation by applying a simple Hubbard- 
Stratonovich transformation. The time-local path integral can 
then easily be solved by propagating an equation of motion for 
the reduced density matrix. The auxiliary quantity introduced 
by the transform appears in the dynamics as an additional clas- 
sical time-varying force. The functional integration over the 
auxiliary quantity can be performed stochastically by inter- 
preting it as a noise force. This yields a new and elegantly 
simple algorithm: Generate an ensemble of Gaussian noise 
trajectories with the appropriate spectrum, propagate the sys- 
tem for each noise trajectory, and average over the ensemble. 
In stark contrast to QMC methods, we find a statistical error 
that is virtually independent of time. 

At a microscopic level, dissipation is caused by the interac- 
tion of a non-thermal system with a vast number of environ- 
mental modes at or close to thermal equilibrium. Although 



cumulative effects of these modes interacting with the system 
under study can have a drastic effect on its dynamics, the cou- 
pling to each individual mode is usually weak. This justifies 
the paradigmatic model employed by Caldeira and Leggett 
Jic| ] to describe dissipative quantum systems - a particle in 
an arbitrary potential interacting linearly with a vast number 
of harmonic excitations, 
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The effect of the environmental modes is fully characterized 

2 

by their spectral density J(ui) — § J2u m"u ~ u>v )' 
which takes on the form J{lo) = r/uj in the Ohmic case, r\ 
being the classical friction constant. 

The time-dependent reduced density matrix is then given 
by a double path integral 



p(q U qf,t) = / V[q] / V[q']e-> 

J q\ J qi 



(So[q]-S [q']) 



F[q,q'}. (2) 



Here So[q] is the action of the undamped quantum sys- 
tem, and its interactions with the environment are incorpo- 
rated into a complex-valued influence functional F[q, q'] = 
exp(-<5>'/h - i&'/h) with 



dt' \ dt"(q(t') 

t J to 

xL'(t'-t")(q(t") 
1 



?'(*')) 
-Q'(t")), 



$'W] = ^? dt'(q(t')-q'(t')) 

L J to 

X(q(t') + q'(t')). 



(3) 



(4) 



L'{t) is the real part of the autocorrelation function of the col- 
lective bath coordinate c u x v , averaged over an ensemble 
of free oscillators. 

This influence functional describes a 'factorized' initial 
preparation, i.e. with the particle constrained to an initial po- 
sition qi for times t < t [|ll]] with the environment fully re- 
laxed. Equilibrium correlation functions may be calculated 
by pushing the preparation back to a sufficiently large nega- 
tive time to < ||] and inserting measurement operators at 
times and t. 
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In the following we outline the numerical solution of these 
equations using a novel technique we call the SLED (Stochas- 
tic Liouville Equation for Dissipation) method. SLED is 
based on the same principles as another method we recently 
introduced - the CSQD (Chromostochastic Quantum Dynam- 
ics) method [O]. Whereas CSQD was developed to treat dis- 
crete systems, SLED solves extended systems using a stochas- 
tic Liouvillian formalism. 

The primary obstacle in trying to translate eqs. (||) - (Q) into 
a recursion relation or an equation of motion for p lies in the 
interaction kernel L'(t). Its range is often large, becoming di- 
vergent for T — > 0. This problem of retarded self-interactions 
mediated by L'(t) can be solved, albeit at the cost of intro- 
ducing an additional path variable. The exponential of the 
non-local action g'] can be decomposed into time-local 
phase factors, 



ex V (-&[q,q'}) = 

2?[£] W[$ exp j-t Jdt'^t')[q{t') - q'{t')) 



(5) 



The distribution function W[£] is real and Gaussian, with 
(£(*)£(*0)w = L'(t - t'), and normalizable through the 
condition q] = 0. Formally, this decomposition is just 
a Hubbard-Stratonovich transformation in a function space 
over the interval [to, t]. Equation (|5|) is also equivalent to the 
construction of an influence functional for a classical colored 
noise source [^], and as such, we will interpret the function 
as a noise trajectory and the measure X>[£]W^[£] as the 
probability measure of an ensemble of noise trajectories. For 
each noise trajectory, the reduced density matrix is then prop- 
agated deterministically according to the equation of motion 



ifip = C p + C v p + £[q, p] 



(6) 



where Co is the Liouvillian of the undamped system, and the 
friction term $" gives rise to an additional operator in the Li- 
ouvillian 



(7) 



The spectrum of the noise is determined by the spectral 
density and the temperature. In the Ohmic case, it is given 
by S(u>) = f]u> coth(hf3ui/2) up to an arbitrary but neces- 
sary cutoff frequency u> c . For computational purposes, it 
is advantageous to split off a 'classical' white noise part 
Sc\(u) — 2r}/h(3 from this spectrum because it can be treated 
by a damping term 



C d 



(8) 



that can be incorporated into the deterministic equation of mo- 
tion in the Liouvillian rather than explicitly in the numerically 



generated noise field [ 1 3 1 . The full equation of motion then 
reads 



ihp = C p + C n p + Cap + £[q, p] 



(9) 



where the noise spectrum of £ is S(u>) — S(u>) — S c \(uj). It 
has been shown previously [|l3|] that using eq. ^ without 
the noise force £ is a limiting case that becomes exact only if 
S(u>) — S c \(ui), i.e., for very high temperatures kT 3> hu> c . 
Comparing the exact dynamics expressed in eq. ([)]) to this 
limiting case, we find a remarkable interpretation of our for- 
malism: When a linear environment is treated classically, one 
generally incurs an error due to the missing quantum fluctua- 
tions. But by explicitly substituting an external noise force for 
the quantum fluctuations one actually regains the exact quan- 
tum dynamics at any temperature ! 

In the following we give examples, solving eq. (^) in the 
position representation using the split-propagator technique. 
The alternating direction implicit method is used for the oper- 
ator Cq + £[q, ■} and an explicit method for the remaining term 
C n + Cd- This choice is appropriate as it preserves the trace 
of p. 

Our first test consists in the computation of the equilib- 
rium correlation function S(t) = ^(q(t)q(0) + q(Q)q(t)} of 
the damped harmonic oscillator, shown in Fig. 1, and its re- 
sponse function \(t) = ^- , 9(0)]}. At relatively high 
temperature, kT = 2hui , this calculation converges with only 
a handful trajectories and reproduces an almost classical re- 
sult, as expected in this temperature regime. The initial value 
5(0) corresponds exactly to the thermal expectation value 
{q 2 ) = kT/niLu 2 . At low temperature, kT = 0.2?uvo, the 
fluctuation amplitude is no longer temperature-dependent, but 
is instead indicative of zero-point fluctuations. (Note that the 
initial amplitude for the damped oscillator is slightly reduced 
from the undamped value (q 2 ) = H/2muJo [|l4|]). In the same 
test, the susceptibility \(t) of the quantum harmonic oscillator 
is determined. We find excellent agreement with the known 
analytic result x(i) = Xcl(t)- The high-temperature calcu- 
lation completes within minutes, while the low-temperature 
calculation takes about 10 3 trajectories for a very small abso- 
lute error < 0.007, or 40 CPU hours on a midrange IBM Rise 
processor using a preliminary, non-optimized code. 

For a more elaborate example, we use the quartic double- 
well potential V(q) — q 4 /4 — q 2 , where natural units for 
time, mass, and energy are implied. In the undamped case, 
this system shows complex dynamics including both tunnel- 
ing and vibrational oscillations. Fig. 2 shows the equilibrium 
dynamics of this system for different values of friction and 
temperature. For high friction, rj = 0.5, the system shows 
a monotonous decay, characteristic of incoherent relaxation 
(solid line). When the friction is lowered to rj = 0.05, a qual- 
itative change to complex, anharmonic oscillatory dynamics 
is observed (dotted line). For the same friction, lowering 
the temperature leads to another qualitative change in the dy- 
namics, resulting in damped oscillations with two manifestly 
distinct frequency scales (dashed line). The low-frequency 
component represents coherent tunneling, with a frequency 
related to the splitting between the two lowest energy eigen- 
states of the double well, and the high-frequency part derives 
from higher excited states. 
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In conclusion, we have developed a novel algorithm for the 
dynamics of extended dissipative quantum systems, which is 
highly efficient. It is also shown to be accurate on first princi- 
ples and by explicit numerical tests. We have presented sim- 
ulation data both for the damped harmonic oscillator and the 
dissipative double-well system which show the practical fea- 
sibility and power of this method. Empirically, we find that 
the statistical error of our results is extremely small, asymp- 
totically approaching a constant value for long times. This 
implies a CPU requirement that scales linearly with the sys- 
tem time t, compared to the exponential rise for conventional 
quantum Monte Carlo methods. The SLED algorithm is also 
easily implemented on parallel computers. The system prop- 
agation for different noise realizations can be performed in 
parallel without any dependencies beyond the accumulation 
of final results. 

We are currently extending the SLED method to study 
dynamics on electronically coupled energy surfaces in 
condensed-matter problems. The problem of greatest interest 
is the dynamics of electron transfer reactions in the inverted 
region. To model a bath with an arbitrary cutoff frequency, 
the SLED method has to be extended to treat systems with 
two electronic states coupled to a solvent mode which is in 
turn coupled to an Ohmic bath JT^]. 
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FIG. 1. Equilibrium Fluctuations of the damped harmonic os- 
cillator at T = 2 (solid line) and T = 0.2 (long dash) with 
h — k = m = 1, ljo = 1, rj — 0.2. Short dashes indicate the 
analytic high-temparature solution [|14|] for T — 2. 
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FIG. 2. Equilibrium fluctuations of a dissipative double-well sys- 
tem for different friction constants. Solid line - dynamics showing 
incoherent decay for r\ = 0.5, T = 1. Dotted line - dynamics show- 
ing predominantly anharmonic vibronic oscillations for r\ = 0.05, 
T — 1. Dashed line - dynamics showing predominantly tunneling 
oscillations for 77 = 0.05, T = 0.5. 
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